home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / src / haeberli / libgutil / vcool.c < prev    next >
C/C++ Source or Header  |  1994-08-01  |  11KB  |  369 lines
  1. /*
  2.  *     vcool - 
  3.  *        Scan an image for pixels with RGB values that will give
  4.  *    "unsafe" values of chrominance signal or composite signal
  5.  *    amplitude when encoded into an NTSC or PAL colour signal.  This
  6.  *    happens for certain high-intensity high-saturation colours that 
  7.  *    are rare in real scenes, but can easily be present in synthetic 
  8.  *    images.
  9.  *
  10.  *     Such pixels can be flagged so the user may then choose other
  11.  *    colours.  Or, the offending pixels can be made "safe"
  12.  *    in a manner that preserves hue.
  13.  *
  14.  *    There are two reasonable ways to make a pixel "safe":
  15.  *    We can reduce its intensity (luminance) while leaving
  16.  *    hue and saturation the same.  Or, we can reduce saturation
  17.  *    while leaving hue and luminance the same.  An argument selects
  18.  *    which strategy to use.
  19.  *
  20.  *    exports:
  21.  *        int vcool(short *red, short *green, short *blue, int n);
  22.  *        void vcoolpixrepair(repair);
  23.  *        void vcoolpixgamma(gam);
  24.  *
  25.  *                     <>
  26.  *
  27.  *     Originally written as "ikNTSC.c" by Alan Wm Paeth,
  28.  *    University of Waterloo, August, 1985
  29.  *
  30.  *    Previous incarnation was hot.c by Dave Martindale,
  31.  *      Imax Systems Corp., December 1990
  32.  *
  33.  *     Updated by Dave Martindale, Imax Systems Corp., December 1990
  34.  * 
  35.  *    Mucked with by David Blythe, Silicon Graphics, December 1991
  36.  *
  37.  *                     <>
  38.  *     Compile-time options.
  39.  *
  40.  * DEFPIXGAMMA
  41.  *     Default gamma of the input and output pixels.  This value is 1.7
  42.  *    for GL applications typically.  This should really be obtained
  43.  *    from the graphics system dynamically.
  44.  *
  45.  * NTSC and PAL
  46.  *     Define either NTSC or PAL as 1 to select the colour system.
  47.  *     Define the other one as zero, or leave it undefined.
  48.  *
  49.  * CHROMA_LIM 
  50.  *    is the limit (in IRE units) of the overall chrominance amplitude; 
  51.  *    it should be 50 or perhaps very slightly higher.
  52.  * 
  53.  * COMPOS_LIM 
  54.  *    is the maximum amplitude (in IRE units) allowed for the composite 
  55.  *    signal.  A value of 100 is the maximum monochrome white, and is 
  56.  *    always safe.  120 is the absolute limit for NTSC broadcasting, 
  57.  *     since the transmitter's carrier goes to zero with 120 IRE input 
  58.  *    signal.  Generally, 110 is a good compromise - it allows somewhat 
  59.  *    brighter colours than 100, while staying safely away from the 
  60.  *    hard limit.
  61.  */
  62. #include "math.h"
  63. #include "stdio.h"
  64. #include "vcool.h"
  65.  
  66. #define    DEFPIXGAMMA    (1.7)        /* this is the default gamma on IRIS */
  67.  
  68. #define    NTSC            (1)
  69. #define    PAL             (0)
  70.  
  71. #define    CHROMA_LIM      (50.0)        /* chroma amplitude limit */
  72. #define    COMPOS_LIM      (110.0)        /* max IRE amplitude */
  73.  
  74. #if NTSC
  75. /*
  76.  * RGB to YIQ encoding matrix.
  77.  */
  78. static float code_matrix[3][3] = {
  79.      0.2989,     0.5866,     0.1144,
  80.      0.5959,    -0.2741,    -0.3218,
  81.      0.2113,    -0.5227,     0.3113,
  82. };
  83.  
  84. #define    PEDESTAL    7.5        /* 7.5 IRE black pedestal */
  85. #define    VIDGAMMA    2.2
  86. #endif /* NTSC */
  87.  
  88. #if PAL
  89. /*
  90.  * RGB to YUV encoding matrix.
  91.  */
  92. static float code_matrix[3][3] = {
  93.      0.2989,     0.5866,     0.1144,
  94.     -0.1473,    -0.2891,     0.4364,
  95.      0.6149,    -0.5145,    -0.1004,
  96. };
  97.  
  98. #define    PEDESTAL    0.0        /* no pedestal in PAL */
  99. #define    VIDGAMMA    2.8
  100. #endif /* PAL */
  101.  
  102. #define SCALE    8192            /* scale factor: do floats with int math */
  103. #define MAXPIX     255            /* white value */
  104.  
  105. /*
  106.  *    if input pixels are in linear space #define CORGAMMA    (VIDGAMMA)
  107.  *    if input pixels are in video space  #define CORGAMMA    (1.0)
  108.  */
  109. #define    CORGAMMA    (VIDGAMMA*(defpixgamma/2.2))    
  110.  
  111. /*
  112.  *  pixel decoding and encoding functions, map integer pixel value
  113.  *  to [0,1] and vice versa.  For now assume a simple linear transformation
  114.  */
  115. #define pix_decode(v)    (((float)(v))/MAXPIX)
  116. #define pix_encode(v)    ((int)((v)*MAXPIX+.5))
  117.  
  118. /*
  119.  * gc: apply the gamma correction specified for this video standard.
  120.  * inv_gc: inverse function of gc.
  121.  *
  122.  * These are generally just a call to powf(), but be careful!
  123.  * Future standards may use more complex functions.
  124.  * (e.g. SMPTE 240M's "electro-optic transfer characteristic").
  125.  */
  126. #define gc(x)        powf((x), (float)1.0 / (float)CORGAMMA)
  127. #define inv_gc(x)    powf((x), (float)CORGAMMA)
  128.  
  129. static int tab[3][3][MAXPIX+1]; /* multiply lookup table */
  130. static float chroma_lim;        /* chroma limit */
  131. static float compos_lim;        /* composite amplitude limit */
  132. static long ichroma_lim2;       /* chroma limit squared (scaled integer) */
  133. static int icompos_lim;        /* composite amplitude limit (scaled integer) */
  134. static float defpixgamma;      /* the default gamma ramp is 1.7 */
  135. static int tabmade = 0;
  136. static int repairmode = 0;
  137.  
  138. /*
  139.  * build_tab: Build multiply lookup table.
  140.  *
  141.  * For each possible pixel value, decode value into floating-point
  142.  * intensity.  Then do gamma correction required by the video
  143.  * standard.  Scale the result by our fixed-point scale factor.
  144.  * Then calculate 9 lookup table entries for this pixel value.
  145.  *
  146.  * We also calculate floating-point and scaled integer versions
  147.  * of our limits here.  This prevents evaluating expressions every pixel
  148.  * when the compiler is too stupid to evaluate constant-valued
  149.  * floating-point expressions at compile time.
  150.  *
  151.  * For convenience, the limits are #defined using IRE units.
  152.  * We must convert them here into the units in which YIQ
  153.  * are measured.  The conversion from IRE to internal units
  154.  * depends on the pedestal level in use, since as Y goes from
  155.  * 0 to 1, the signal goes from the pedestal level to 100 IRE.
  156.  * Chroma is always scaled to remain consistent with Y.
  157.  */
  158. static build_tab()
  159. {
  160.     float    f;
  161.     int pv;
  162.  
  163.     for (pv = 0; pv <= MAXPIX; pv++) {
  164.         f = SCALE * gc(pix_decode(pv));
  165.         tab[0][0][pv] = (int)(f * code_matrix[0][0] + 0.5);
  166.         tab[0][1][pv] = (int)(f * code_matrix[0][1] + 0.5);
  167.         tab[0][2][pv] = (int)(f * code_matrix[0][2] + 0.5);
  168.         tab[1][0][pv] = (int)(f * code_matrix[1][0] + 0.5);
  169.         tab[1][1][pv] = (int)(f * code_matrix[1][1] + 0.5);
  170.         tab[1][2][pv] = (int)(f * code_matrix[1][2] + 0.5);
  171.         tab[2][0][pv] = (int)(f * code_matrix[2][0] + 0.5);
  172.         tab[2][1][pv] = (int)(f * code_matrix[2][1] + 0.5);
  173.         tab[2][2][pv] = (int)(f * code_matrix[2][2] + 0.5);
  174.     }
  175.  
  176.     chroma_lim = (float)CHROMA_LIM / (100.0 - PEDESTAL);
  177.     compos_lim = ((float)COMPOS_LIM - PEDESTAL) / (100.0 - PEDESTAL);
  178.  
  179.     ichroma_lim2 = (int)(chroma_lim * SCALE + 0.5);
  180.     ichroma_lim2 *= ichroma_lim2;
  181.     icompos_lim = (int)(compos_lim * SCALE + 0.5);
  182. }
  183.  
  184. void vcoolpixgamma(gam)
  185. float gam;
  186. {
  187.     defpixgamma = gam;
  188.     build_tab();
  189.     tabmade++;
  190. }
  191.  
  192. void vcoolrepairmode(repair)
  193. int repair;
  194. {
  195.     repairmode = repair;
  196. }
  197.  
  198. int vcool(short *red, short *green, short *blue, int n)
  199. {
  200.     int    r, g, b;
  201.     int    y, i, q;
  202.     long    y2, c2;
  203.     float        pr, pg, pb;
  204.     float        py;
  205.     float    fy, fc, t, scale;
  206.     static int    prev_r = 0, prev_g = 0, prev_b = 0;
  207.     static int    new_r, new_g, new_b;
  208.     int nhot = 0;
  209.  
  210.     if (!tabmade)
  211.         vcoolpixgamma(DEFPIXGAMMA);
  212.     if (!repairmode) 
  213.         vcoolrepairmode(REDUCE_SAT);
  214.  
  215.     while(n--) {
  216.         r = *red++;
  217.         g = *green++;
  218.         b = *blue++;
  219.  
  220.         /*
  221.          * Pixel decoding, gamma correction, and matrix multiplication
  222.          * all done by lookup table.
  223.          *
  224.          * "i" and "q" are the two chrominance components;
  225.          * they are I and Q for NTSC.
  226.          * For PAL, "i" is U (scaled B-Y) and "q" is V (scaled R-Y).
  227.          * Since we only care about the length of the chroma vector,
  228.          * not its angle, we don't care which is which.
  229.          */
  230.         y = tab[0][0][r] + tab[0][1][g] + tab[0][2][b];
  231.         i = tab[1][0][r] + tab[1][1][g] + tab[1][2][b];
  232.         q = tab[2][0][r] + tab[2][1][g] + tab[2][2][b];
  233.  
  234.         /*
  235.          * Check to see if the chrominance vector is too long or the
  236.          * composite waveform amplitude is too large.
  237.          *
  238.          * Chrominance is too large if
  239.          *
  240.          *    sqrt(i^2, q^2)  >  chroma_lim.
  241.          *
  242.          * The composite signal amplitude is too large if
  243.          *
  244.          *    y + sqrt(i^2, q^2)  >  compos_lim.
  245.          *
  246.          * We avoid doing the sqrt by checking
  247.          *
  248.          *    i^2 + q^2  >  chroma_lim^2
  249.          * and
  250.          *    y + sqrt(i^2 + q^2)  >  compos_lim
  251.          *    sqrt(i^2 + q^2)  >  compos_lim - y
  252.          *    i^2 + q^2  >  (compos_lim - y)^2
  253.          *
  254.          */
  255.         c2 = (long)i * i + (long)q * q;
  256.         y2 = (long)icompos_lim - y;
  257.         y2 *= y2;
  258.         if (c2 <= ichroma_lim2 && c2 <= y2)     /* no problems */
  259.             continue;
  260.         /*
  261.          * Pixel is hot, choose desired strategy
  262.          */
  263.         if (repairmode == FLAG_HOT) {
  264.         /*
  265.          * Set the hot pixel to black to identify it.
  266.          */
  267.         r = g = b = 0;
  268.         } else {
  269.         /*
  270.          * Optimization: cache the last-computed hot pixel.
  271.          */
  272.         if (r == prev_r && g == prev_g && b == prev_b) {
  273.             r = new_r;
  274.             g = new_g;
  275.             b = new_b;
  276.             goto update;
  277.         }
  278.         prev_r = r;
  279.         prev_g = g;
  280.         prev_b = b;
  281.  
  282.         /*
  283.          * Get Y and chroma amplitudes in floating point.
  284.          *
  285.          * If your C library doesn't have fhypot(), just use
  286.          * fhypot(a,b) = sqrt(a*a + b*b);
  287.          *
  288.          * Then extract linear (un-gamma-corrected) floating-point
  289.          * pixel RGB values.
  290.          */
  291.         fy = (float)y / SCALE;
  292.         fc = fhypot((float)i / (float)SCALE, (float)q / (float)SCALE);
  293.  
  294.         pr = pix_decode(r);
  295.         pg = pix_decode(g);
  296.         pb = pix_decode(b);
  297.  
  298.         /*
  299.          * Reducing overall pixel intensity by scaling
  300.          * R, G, and B reduces Y, I, and Q by the same factor.
  301.          * This changes luminance but not saturation, since saturation
  302.          * is determined by the chroma/luminance ratio.
  303.          *
  304.          * On the other hand, by linearly interpolating between the
  305.          * original pixel value and a grey pixel with the same
  306.          * luminance (R=G=B=Y), we change saturation without
  307.          * affecting luminance.
  308.          */
  309.         if (repairmode != REDUCE_SAT) {
  310.             /*
  311.              * Calculate a scale factor that will bring the pixel
  312.              * within both chroma and composite limits, if we scale
  313.              * luminance and chroma simultaneously.
  314.              *
  315.              * The calculated chrominance reduction applies to the
  316.              * gamma-corrected RGB values that are the input to
  317.              * the RGB-to-YIQ operation.  Multiplying the
  318.              * original un-gamma-corrected pixel values by
  319.              * the scaling factor raised to the "gamma" power
  320.              * is equivalent, and avoids calling gc() and inv_gc()
  321.              * three times each.
  322.              */
  323.             scale = chroma_lim / fc;
  324.             t = compos_lim / (fy + fc);
  325.             if (t < scale)
  326.                 scale = t;
  327.             scale = powf(scale, (float)CORGAMMA);
  328.  
  329.             r = pix_encode(scale * pr);
  330.             g = pix_encode(scale * pg);
  331.             b = pix_encode(scale * pb);
  332.         } else {
  333.             /*
  334.              * Calculate a scale factor that will bring the pixel
  335.              * within both chroma and composite limits, if we scale
  336.              * chroma while leaving luminance unchanged.
  337.              *
  338.              * We have to interpolate gamma-corrected RGB values,
  339.              * so we must convert from linear to gamma-corrected
  340.              * before interpolation and then back to linear afterwards.
  341.              */
  342.             scale = chroma_lim / fc;
  343.             t = (compos_lim - fy) / fc;
  344.             if (t < scale)
  345.                 scale = t;
  346.  
  347.             pr = gc(pr);
  348.             pg = gc(pg);
  349.             pb = gc(pb);
  350.             py = pr * code_matrix[0][0] + pg * code_matrix[0][1]
  351.                 + pb * code_matrix[0][2];
  352.             r = pix_encode(inv_gc(py + scale * (pr - py)));
  353.             g = pix_encode(inv_gc(py + scale * (pg - py)));
  354.             b = pix_encode(inv_gc(py + scale * (pb - py)));
  355.         }
  356.  
  357.         new_r = r;
  358.         new_g = g;
  359.         new_b = b;
  360.         }
  361. update:
  362.         red[-1] = r;
  363.         green[-1] = g;
  364.         blue[-1] = b;
  365.         nhot++;
  366.     }
  367.     return nhot;
  368. }
  369.